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Abstract 

We explain how, given a plane algebraic curve C : f{x, y) = 0, xi G C 
not a singularity of y w.r.t. x, and e > 0, we can compute 5 > 0 such that 
\yj{xi) — yj{x2)\ < s for all holomorphic functions yj{x) which satisfy 
f(x,yj{x)) = 0 in a neighbourhood of x\ and for all X 2 with \xi — X 2 I <5. 
Consequently, we obtain an algorithm for reliable homotopy continuation 
of plane algebraic curves. As an example application, we study continuous 
deformation of closed discrete Darboux transforms. 

Moreover, we discuss a scheme for reliable homotopy continuation of 
triangular polynomial systems. A general implementation has remained 
elusive so far. However, the epsilon-delta bound enables us to handle the 
special case of systems of plane algebraic curves. The bound helps us 
to determine a feasible step size and paths, which are equivalent w.r.t. 
analytic continuation to the actual paths of the variables but along which 
we can proceed more easily. 


1 Motivation 


In many geometric problems, variables depend analytically on some parameter. 

If we want to analyze and experiment with these problems using interactive soft¬ 
ware, whenever the user continuously modifies the parameter, we must update 
the dependent variables accordingly. For many applications, in doing so, the 
analytical relationship between variables and parameter should be preserved at 
all times. Therefore we need reliable algorithms for analytic continuation. 

for example the following problem of discrete differential geometry (IJloff- 
Section 2.6). Let there be a regular discrete curve 7 in CP^, i.e. 
chain with distinct vertices 7 o, 7 i,..., 7 n S CP^. We define the 
discrete Darboux transform 7 of 7 with initial point 70 S CP^ and parameter 
fj, G C as follows: for all j = 1,2,... ,n, let y G CP^ be the unique point for 
which the cross-ratio 
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It can be shown that 7 j-i is mapped to 7 ^ by a unique Mobius transformation, 
which depends only on 7 j-i, 7 j, and js, but not on 7 _,-. Hence, there exists a 
unique Mobius transformation M depending on 70 , 71 , ..., 7 n, and /i, which 
maps an initial point 70 to the corresponding last point 7 „ of 7 . Consequently, 
for every choice of /r G C, there are two choices of initial point 70 (counted 
with multiplicity) such that 7 is a closed polygonal chain. These are exactly the 
fixed points of M or, in other words, the roots of the characteristic polynomial 
of M. The vanishing of the characteristic polynomial establishes an algebraic 
(particularly analytical) relationship between fj, and 70 . 

If we want to study closed Darboux transforms of a discrete curve 7 for vary¬ 
ing parameter ^ using interactive software, then we must analytically continue 
70 - Otherwise we may observe sudden jumps of 70 under continuous movement 
of fi, which have no mathematical justification. 

In practice, of course, we cannot modify a parameter continuously. Instead, 
we obtain a series of parameter values at a series of discrete points in time. 
We do not know how the parameter moves between sample points. A natural 
approach would be to interpolate linearly between consecutive parameter values 
(using a time parameter in the unit interval). However, the segment between 
parameter values may contain singularities beyond which analytic continuation 
becomes impossible. Thus it seems reasonable to analytically continue along the 
polygonal chain of parameter values as long as this is possible, and to deviate 
from that path otherwise. Such a deviation can still be interpreted as a linear 
interpolation between consecutive parameter values if we let the time parameter 
run from 0 to I on an arbitrary path through the complex plane instead of 
restricting it to the unit interval. 

This is the paradigm of ‘complex detours’ invented by Kortenkamp and 


Richter-Gebert for their interactive geometry software Cinderella (Kortenkamp 


and Richter-Gebert, 20061. It is described in more detail in (Kortenkamp, 1999 


esp. Chapter 7; ?KortenkampRichterGebert2001b Kortenkamp and Richter- 


Gebert, 2002 1 . Essentially the same concept was conceived in the context 


of homotopy continuation by Morgan and Sommese (Morgan and Sommese, 


19871, who later named it the ‘gamma trick’ (Sommese and Wampler, 2005 


Lemma 7.1.3 on p. 94). 

Once we have chosen a path for the parameter, we must determine the right 
value of the dependent variable at consecutive sample points. How this can be 
achieved may in fact be relatively easy to see for us—^just determine values in a 
way such that there are no jumps—but hard to see for an algorithm. The tracing 
problem of dynamic geometry, i.e. tracing the positions of dependent elements 
of a geometric construction under movement of a free element, is NP-complete 
already for constructions that only involve points, lines through two points, in¬ 
tersection of lines, and angle bisectors (Kortenkamp and Richter-Gebert, 2002 1 . 

The interactive geometry software Cinderella currently uses a heuristic for 
path following. Most homotopy continuation methods use a predictor-corrector 
approach, which is generally also heuristic. For an overview of homotopy 
continuation methods, consider the books by Allgower and Georg (1990) or 
Sommese and Wampler (20051. Lately, certified homotopy continuation meth¬ 


ods have emerged (Beltran and Leykin, 2012 2013 Hauenstein and Sottile 


2012 


1986). 


Hauenstein et ah, 2014). They are based on Smale’s alpha theory (Smale, 


In what follows, we derive a certified algorithm for analytic continuation 
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of plane algebraic curves based on the following simple observation: Due to 
continuity, if the parameter moves little, so does the dependent variable. Hence, 
if we take small enough steps along the parameter path, we can choose the right 
value of the dependent variable based on proximity. As an application, we return 
to the example of continuous deformation of closed discrete Darboux transforms. 
Moreover, we show how the algorithm generalizes to systems of plane algebraic 
curves. A comparison with other approaches demonstrates the practicability of 
our algorithms. 

2 Computing an epsilon-delta bound for plane 
algebraic curves 

Theorem 2.1. LetC: f{x,y) = 0 be a complex plane algebraic curve, where 

n 

f{x,y) = ^ak{x)y'^-'^ 

k=0 

is a polynomial of degree n in y whose coefficients ak{x) are polynomials in 
X. Let Xi € C be a point in the complex plane at which neither the leading 
coefficient ao(x) nor the discriminant of f(x,y) w.r.t. y vanish. Then for every 
e > 0, we can algorithmically compute <5 > 0 .such that 

\yj{xi) - yj{x 2 )\ < e 

for all holomorphic functions yj{x), j = 1,2,... ,n, that satisfy f(x, yj(x)) = 0 
in a neighbourhood of X\ and for all X 2 with \xi — X 2 \ < 5. 

Remark 2.2. How does [Theorem 2TT] help us to perform analytic continuation? 
Let e be half the minimal distance between the y-values at xi. Then for any X 2 
less than 5 away from xi the following holds: The j/-value yj{x 2 ), which results 
from analytic continuation of yj{x) along the segment from xi to X 2 , is closer 
to yj{xi) than to any other y-value at Xi. In other words, 5 provides an upper 
bound for the step width of parameter x such that we may match y-values on 
the same branch based on proximity. 

Our plan for the proof of |Theorem 2.1| is as follows: We will see that there is an 
upper bound of <5 depending on 

1. the radius of convergence of the Taylor expansion of yj{x) at Xi, 

2 . the modulus of the derivative of yj{x) at Xi, 

3. the maximum modulus of yj{x) on a circle centred at Xi, 

for j = 1, 2,..., n, respectively. We derive a formula for that upper bound and 
then compute bounds for its ingredients. To this end, we need the following 
lemmas. 

Lemma 2.3. Let U G C be an open subset of the complex plane, and let 

j/j :[/—>■ C 
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he holomorphic. Taylor expansion of yj around xi G U yields 

yj{x2) = Vjixi) + ix2 - Xi)yj{Xi) + ix2 - XifR{x2), 

for all X 2 G C such that \x 2 — xi\ < p and sufficiently small p > 0. The 
remainder R{x 2 ) satisfies 


\R{x2)\ < 


M 

p{p- \X2 - a;i|) 


where 


M= max lyJxi + pe^*)\. 

[0,27r] 


Lemma 2.3 is a standard result of complex analysis (Ahlfors, 1979, p. 124-126), 


which we therefore do not prove here. 


Lemma 2.4 (implicit differentiation). Let f{x,y) he a complex polynomial. Let 
U C C he an open subset of the complex plane. Let pj : U ^ C he a holomorphic 
function that satisfies f{x,yj(x)) = 0 for all x G U. Then for all xi G U with 
fy{xi,yj{xi)) 0 it follows that 


y'jixi) = 


f^{xi,yj{xi)) 

fy{xi,yj{xi))' 


Proof. By the chain rule, the total differential of f{x,yj{x)) = 0 w.r.t. x is 


Df{x,yj{x)) = fo,(x,yj{x)) + fy{x,yj{x)) ■ yl{x) = 0. 


Therefore 


y'j{xi) 


f^{xi,yj{xi)) 

fy{xi,yj{xi))' 


□ 


Lemma 2.5 (Fujiwara (1916 Inequality 3 on p. 168)). Consider a polynomial 


P{x) = ^ ttkX^ 
k=0 


of degree n with complex coefficients Ok G C, k = 0,1,... ,n. Then all x G C 
with p{x) = 0 satisfy 


|a;| < 2max 



: fc = 1, 



Proof. Consider the inequality 


\p{x)\ > |ao||a;r-^|afc||a;r . 


( 1 ) 


The RHS of Q is positive if 

loolla^r > 2 '"|afc||a;|"“'', k = 1,2,...,n, 


4 










because then 


laoll^r > (1 - 2-")|ao||a:r = ^ > E 

Hence, |p(a:)| > 0 if 


n—k 


fc=l 


k^l 


\x\ > max i 2 


and thus 


Ixl < 2 max ■ 


: fc = 1 ,..., n 


for all zeros a; e C of p(x). 


□ 


Lemma 2.6 (bounds for trigonometric polynomials). Consider a trigonometric 
polynomial of degree n of the form 


p{xi + pe“) = ^ afc(a;i + pe“) 


i—k 




Then 


|p(a;i +pe“)| < ^ |afe|(|a;i| + |p|) 


n—k 




Moreover, if the zeros Xi, X 2 , ■ ■ ■, Xn of p(x) satisfy \xk — Xi\ > p then 


\p[xi + / 9 e“)| > |ao| ]^(|a;fe - xi\ - p) > 0 . 
k=0 

Proof. The upper bound follows from the triangle inequality. The lower bound 
follows from the factorization 

n 

p{xi + pe‘*) = ao n (^1 + ~ 

k=0 


and the fact that \xi + pe‘‘ — Xk\ > \p — \xk — xi\\. Note that the lower bound 
is positive by the assumptions that \xk — xi\ > p and that p has degree n, i.e. 

Uq 7^ 0. □ 


Proof of \ Theorem 2. 1\ Let yj{x), j = denote the holomorphic func¬ 

tions that satisfy f{x, yj{x)) = 0 in a neighbourhood of xi. By Lemma 2.3 

yj{x2) = Vjixi) P (X2 - Xi)y'Axi) P (X2 - Xif Rj{X2) ( 2 ) 


for all X 2 S C such that |x 2 — Xi| < p and sufficiently small p > 0. If we bring 
Pj(xi) to the LHS of take the absolute value on both sides, and apply the 
triangle inequality, we see that 


|j/j(a^l) - yj{X2)\ = \X2 - XiWy'jixi) -f (X2 - Xi)Rjix2)\ 

< \X2 - Xi|(|p'(xi)| -I- |X2 - Xi||i?j(x2)|) 

= |i?j(X 2 )||x 2 - XiP -k |p'(xi)||X 2 - Xi|. ( 3 ) 
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Hence, under the above assumptions, 

\Rj{x2)\\x 2 - Xi\^ + \yj{xi)\\x2 - a;i| - e < 0 . ( 4 ) 

is a sufficient condition for \yj{xi) — yj{x 2 )\ < £• 

The LHS of Q is strictly increasing in |y'(xi)| and \Rj{x 2 )\- Therefore, if 
we plug in the bounds 

\yjixi)\ < max|y'(a:i)| =: Y (5) 

j 3 


and 


\Rj{X2)\ < 


M 


p{p- |a:2 -a;i|) 

(see Lemma 2.31 into we obtain a stronger sufficient condition for 

\yj{xi) - yj{x2)\ <£, 


namely 


-7- ^ -rTk2 -Xi\^+ Y\x2 - 

pip - \X2-Xi\) 

M\x2 - XiP + p{p - \X2 - Xi\){Y\x2 
(M - pT)|a;2 - + p{pY + £)|a;2 - 


Xil — £ < 0 

— xi\ — e) <0 
- xi\ — ep^ < 0. 


( 6 ) 


How we can transform (§ into a sufficient bound on \x 2 — xi\ depends on the 
sign of M — pY. 

First case: M — pY > 0. The LHS of (|^ describes a smile parabola in 
\x 2 — xi\ with a positive and a negative root. Since \x 2 — xi\ > 0, we need only 
bound \x 2 — xi\ from above by the positive root, i.e. 


-p{pY + £) + J p'^ipY + £)^ + 4(M - pY)ep^ 

1 "^ - "‘I - 

p - {pY + 

^ 2(M - pY) ■ 

Second case: M — pY < 0. The LHS of (§ describes a frown parabola in 
\x 2 — xi\ with one root greater than p and one root between 0 and p. Since 
\x 2 — xi\ < p by definition, we need only bound \x 2 — xi\ from above by the 
smaller root, i.e. 


\X2 -Xi\< 


p{pY + £) - P^J(pY + e)^ - 4:{pY - M) 
~ 2{pY - M) 

— £)^ + 4£M — {pY + e) 

2{M - pY) 


Third case: M — pY = 0. The LHS of ^ reduces to 

p(pF + £)|a;2 - a;i| - £p^ < 0 \x2-xi\< 


ep 


pY Ye 
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This bound is asymptotically equivalent to the previous bounds as M approaches 
pY. Altogether, we thus arrive at the sufficient bound 


pi\J {pY — e)^ + 4eM — {pY + e 

- 2(M - fY) - 

The RHS of Q has the expected qualitative behaviour: It is strictly increasing 
in e and p, and strictly decreasing in M and Y. 

It remains to be shown that we can compute bounds for the ingredients p, 
Y, and M of ([^ . 

Lemma 2.3| (and thus our argument) is valid if and only if p is smaller than 
the radius of convergence of the Taylor expansion of yj{x). Therefore, we must 
choose p smaller than the distance between xi and the singularities of yj{x), 
j = 1, 2,..., n. Recall that yj{x) satisfies f{x, yj{x)) = 0 in a neighbourhood of 
xi, where 

n 

f{x,y) = '^ak{x)y'^-^. 
k=0 

In particular, p must be smaller than the distance between Xi and the zeros of 
ao(x). The zeros of ao(x) are exactly the poles of yj(x). The remaining finite 
singularities of yj{x) are exactly the finite ramification points of yj{x). These 
are zeros of the discriminant of f{x,y) w.r.t. y. Hence, we may choose any 



p < min{|a;i - a;|: qq^x) ■ Ay{f{x,y)){x) = 0}, 

where Ay{f{x,y)){x) denotes the discriminant of / w.r.t. y. 
We can compute 


Y = max |y'(xi)| = max 


fx{xi,yj{xi)) 


fy{xi,yj{xi)) 


by [Lemma 2.4 Note that the denominator does not vanish by the assumption 
that xi is not a zero of the discriminant of f{x,y) w.r.t. y. 

Therefore, M remains to be computed or bounded from above. To that end, 
we can apply [Lemma 2.'5]to 


f{x,yj{x)) =J2°'k{x)yj{x)^, 


fc =0 


interpreted as a polynomial in yj{x). By our choice of p, the leading coefficient 
ao{x) does not vanish for all x with |x — xi| < p. For those x and for all 
j = 1, 2,..., 71, [Lemma 2.5[yields 


\yj{x)\ < 2max 


ak{x) 

ao{x) 



Consequently, 


M <2 max 

tG[ 0 , 27 r] 


ak{xi + pe**) 
ao(xi + pe**) 
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By 


Lemma 2.6 


we can compute upper bounds hk of maX(g[o, 27 r] \o-k{xi + pe'*)| 


and a lower bound oq > 0 of mintg[o^ 27 r] |ao(2;i + pe'*)|, which are much easier 
to compute than these extreme values. 

The zeros of ao(a:) and of Ay{f{x,y))(x) can be computed (at least to ar¬ 
bitrary precision) using a root-finding algorithm. Similarly, the values yj(xi), 
j = 1, 2 ,..., n, can be computed (at least to arbitrary precision) by solving 


f{xi,yj{xi)) = 0 


for yjixi). 

Let us summarize our argument: We may choose 


d = 


p {pY - e)^ -I- AeM - {pY + e)^ 


2(M - pY) 


( 8 ) 


where 


p < min{|a;i — a;|: ao(x 

f^{xi,yj{xi)) 


Y := max 
j 


fy{xi,yj{xi)) 


^y{fi.x,y)){x) = 0}, 


M := 2 max ( ^ 

k \ Oq 


□ 


Remark 2.7. For Theorem 2.1| to hold, f{x,y) needs neither be irreducible nor 
square-free. However, if f{x,y) is not square-free, the discriminant may vanish 
identically and the epsilon-delta bound is no longer useful. If f{x,y) is square- 
free but not irreducible, the epsilon-delta bound for y-values on one irreducible 
component may be smaller than necessary due to the influence of zeros of the 
discriminant of other irreducible components. 


3 Certified homotopy continuation of plane al¬ 
gebraic curves 

Theorem 2.1 [ enables us to solve the following problem: 

Problem 3.1. Consider a plane algebraic curve 

C- fix,y) = 0 . 


Let x: [0,1] —>■ C, t !->■ x{t) be a monotonic (distance non-decreasing) path, i.e. 
|a;(0) — x{ti)\ < |x(0) — x{t 2 )\ for 0 < < ^2 < 1- 


Let j/(0) S C satisfy f{x{0),y{0)) = 0. If analytic continuation of y along x{t) 
is possible, determine the value ?/(l) that results from initial value y(0) under 
analytic continuation of y along x(t). 


The algorithm for Problem 3.1| follows from [Remark 2.2 


Algorithm 3.2. Let f{x,y), x(t), and y{0) be defined as in Problem 3.1 


1. Let T = 0. 





















2. While T < 1, 

(a) Let e be half the minimum distance between the y with 

fixiT),y)=0. 


(b) Compute S by the epsilon-delta bound of Theorem 2.1 


(c) Use bisection to maximize T* G [T, 1] such that \x{T) — x{T*)\ < 5. 

(d) Let y{T*) be the y with f{x{T*),y) = 0 closest to y{T). 

(e) Let T = r*. 

3. Output 2 /( 1 ) and stop. 


Case study: continuous deformation of closed 
discrete Darboux transforms 


[Algorithm 3. 2| shows how the epsilon-delta bound can be used for certified ho- 
motopy continuation of plane algebraic curves. In this section, as an example 
application, let us return to the closed discrete Darboux transform introduced 
in [section 11 


We generally follow the exposition of Hoffmann (2009 Section 2.6) but use 


a slightly different definition of cross-ratio. (A value y of our cross-ratio corres¬ 
ponds to a value 1 — /r of the cross-ratio in (Hoffmann, 2009 Section 2.6) and 
vice versa.) 

Recall the definition of discrete Darboux transform: 

Definition 4.1 (discrete Darboux transform). Let 7 be a regular discrete curve 
in CP^ with vertices 70 , 71 , ..., 7 n G CP^. We choose an initial point 70 G CP^ 
and prescribe a cross-ratio /r G C. The discrete Darboux transform of 7 with 
initial point % and parameter p, is the unique discrete curve 7 whose vertices 
Ip i = 1 , 2 , ...,n, satisfy 


(7j-i,7i;7i,7i-i) := 


(7j-i - 7j)(7j- - U-i) 

(7j-i - 7i-i)(7j - 1]) 


= A*- 


Lemma 4.2. Let a,b,d G CP^ be in general position. For every p G C, there 
exists a Mobius transformation depending on a, b, and p that maps d to c G CP^ 
such that (a, &; c, d) = p. 

Proof. Consider the Mobius transformation 


M: X 


x — a 

X — 6 ’ 


which maps a, b, and d to 0, 00 , and d' respectively. The cross-ratio is invariant 
under Mobius transformations. Hence, if we denote the image of c under M as 
d, we want that 


( 0 , 00 ; d, d') = 


(0 — c')(oo — d') d 
(0 — (i')(oo — d) d' ^ 
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We define the Mobius transformations 


N: d! ^ d = iidf, 


M-^ 


X' !->■ 


bx' — a 
x' — 1 ' 


Then the Mobius transformation 


M ^ o N o M: 


(fib — a)d — (fi — l)ab 
(fi — l)d + b — fia 


maps d G CP^ to c G CP^ such that (a, b; c, d) = fi. □ 

Note that (M~^ o N o M)(a) = a and (M~^ o N o M)(b) = 6, independent of fi. 


Lemma 4.3. There exists a Mobius transformation depending on 70, 71, 

7„, and fi that maps an initial point % of a discrete Darboux transform of ^ 
with parameter fi to the corresponding end point 7 „. 


Proof. By Lemma 4.2 there exist Mobius transformations Mj, j = 1,2,... ,n, 
depending on 7^-1, 'jj, and fi that map 7^-1 to fj. Therefore their composition 
Mn o Mn-1 o • • • o Ml is a Mobius transformation depending on 70, 71, ..., 7„, 
and fi that maps % to 7„. □ 


Remark 4 . 4 . A discrete Darboux transform 7 is closed if and only if its initial 
point 7o is a fixed point of the Mobius transformation of [Lemma 4 . 3 | The 
Mobius transformation of ILemma 4.31 is of the form 

ax + h 


where a, b, c, and d are polynomials in fi with complex coefficients depending 
on 7o, 71, ..., 7n. Its fixed points are the roots of the equation 


(cx + d)x — (ax + b) = cx^ + (d— a)x — 6 = 0. 


This equation is quadratic in x. Its degree in fi increases with the number of 
points of 7. Equivalently, in homogeneous coordinates, the fixed points are the 
eigenvectors of matrix 

a b\ 
c d) ■ 

Example 4.5. As a simple but interesting enough example, consider the closed 
discrete curve 7 spanned by the fifth roots of unity, 

J= 0 ,I,..., 5 . 


The relationship between fi and the initial point 70 of a closed discrete Darboux 
transform 7 of 7 is governed by the equation 


((-3 + fi^ + 6 fi- 3 - 75) 7^ + ((-2 - 4^5) Ai + I + V5) 70 
+ (-3 + 75) + fi/r - 3 - 75 (1 - Ai) = 0. 


( 9 ) 


Equation ([^ is quadratic in 70, cubic in fi, and has total degree 5 . For almost 
every value of fi, exactly two values of 70 satisfy the equation. The only excep¬ 
tions are fi = 1, where all values of 70 satisfy the equation, and fi = and 
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Figure 4 . 6 : continuous deformation of a closed discrete Darboux transform 


fj, = oo, which are ramihcation points of 70, i.e. points where there is only one 
value of 7o, which is a root of multiplicity 2 of (©• 

The discrete Darboux transform 7 of 7 with initial point 70 = 71 and para¬ 
meter /r = 0 is identical to 7 up to a rotation by 27 r/ 5 , i.e. 

7j-i = • 7j-i = 7j 


for all j = 0 , 1 ,..., n. Particularly, the discrete curve 7 is closed. 

We would like to examine how the closed Darboux transform 7 behaves when 
/r makes two full anticlockwise turns around the ramification point on a 

circle through the origin centred at ^ ^ /2 + ^ 


1000 ' 


[Figure 4 T 6 | attempts to illustrate the behaviour of the closed Darboux trans¬ 
form 7 under the aforementioned motion of fj,. (It is notoriously difficult to 
reproduce dynamic behaviour on paper. A video of the experiment is available 
in the supplementary material for this article.) We can read Figure T^ in two 
different ways: 

Firstly, the left image shows the movement of 70 (grey points) as fjL (white 
points) completes one full circle. Then the right image shows the movement of 
7o (grey points) as /i (white points) completes another full circle. The position 
of the Darboux transform (black) after one turn of /r is identical to the initial 
position up to rotation. The hnal position of the Darboux transform after the 
second turn is absolutely identical to the initial position. 

Secondly, the left image shows the movement of one choice of 70 such that 
7 is closed as /r completes one full circle. The right image shows how the other 
choice of 70 such that 7 is closed moves at the same time. After one turn of /i 
we reach the initial position up to interchanged choices of 70. (In the left image, 
% moves from 71 to 74 while in the right image, 70 moves from 74 to 74.) After 
another turn of /i, the two choices of 70 reach their initial positions again. 

Following [Remark 2.2 the steps of /r are chosen according to the epsilon- 
delta bound of Theorem 2.1 for (|^ with e half the distance between the two 
choices of 70. As we expect, the closer /i approaches the singularity the smaller 
the steps become. At its rightmost position, fj, is only away from the 
singularity. It takes 127 steps until fj, completes one full circle. 
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Remark 4 . 7 . If we want to prevent jumps, /i and 70 cannot both be freely 
movable, i.e. we cannot let /i and 70 interchange their roles as movable and 
dependent point. Otherwise, we can force a jump as follows: We move parameter 
p, to /r = 1. At the same time, according to ([^, the two possible initial points 
of 7 move to the origin and the point at infinity, respectively. Without loss of 
generality, we assume that 70 moved to the origin. Note that /i = 1 describes an 
irreducible component of the plane algebraic curve This means that if we 
remove 70 from the origin, /r will simply rest at /r = 1. Then we cannot move /i 
without a jump of 70 because in order to move continuously, 70 would initially 
have to be arbitrarily close to the origin (or the point at infinity). 


Remark 4 . 8 . Floating point arithmetic introduces rounding errors into the com¬ 
putation of 7o. This has a peculiar effect: If 7 is closed, we know that 70 must 
be one of the two fixed points of the Mobius transformation M that maps the 
initial point % of 7 to its last point 7^. In general, a Mobius transformation has 
one attracting and one repelling fixed point. When the fixed point is repelling, 
any numerical error in its position is amplified by Mobius transformation M. 
Therefore, the closed Darboux transform may (numerically) no longer be closed 
when computed naively. We have observed (see Figure that we can move 
from one choice of % to the other by moving /r around a ramification point. The 
natural domain for the map c: /x 1—> 70 is a Riemann surface. Different choices of 
7o correspond to different branches of the Riemann surface. When we compute 
the vertices of 7, we step by step compute M o c = Mn o M^-i o • • • o Mi o c. 
Note that function M o c is an example of a function on a Riemann surface that 
is numerically stable on one branch and numerically unstable on the other. 

Luckily, we can stabilize the computation by considering the inverse Mobius 
transformation M~^. A repelling fixed point of a Mobius transformation is an 
attracting fixed point of its inverse. We can step by step compute M~^ o c = 


o M2^ o ■ ■ ■ o M~^ o c to obtain the vertices of 7 in reversed order. Since 
the cross-ratio of A, B, C, D satisfies (A, B] C, D) = {B, A; £>, C), we need only 
change our algorithm very little in order to obtain M~^ = Mi o M2 o ■ ■ ■ o M„ 
instead of M = Mn o Mn-i o • • • o Mi; we only need to reverse the order of the 
vertices of 7 before we compute the Mobius transformations. 

Besides, we can determine whether 70 approximates an attracting fixed point 
of M by considering the derivative of M at 79. The fixed point near 70 is 
attracting if the absolute value of the derivative is smaller than 1. 


5 Towards homotopy continuation of triangular 
systems of polynomials 


In this section, we discuss a scheme for certified homotopy continuation of trian¬ 
gular systems of polynomials. A general implementation has remained elusive 
so far. However, we later follow the same scheme when we derive an algorithm 


for certified homotopy continuation of systems of plane algebraic curves (Al¬ 


gorithm 6 . 4 ). 


Problem 5.1. Consider a triangular system of polynomials, without loss of 
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generality 

Pi(a;o,a;i) = 0, 

P2ixo,Xi,X2) = 0 , 

( 10 ) 

Pnixo,Xi, ...,Xn)=0- 

Let a;o(0), a;i(0), a;„(0) be initial values that satisfy the system of equa¬ 

tions and let a;o(l) be a target value for variable xq, i.e. a value to which xq 
should move continuously. We define function xo(t) as a parameterization of 
the segment between Xo(0) and 2:9 (1), 


xo(t) = (1 - t)a;o(0) -f teo(l). 


By analytic continuation w.r.t. t S [0,1] we can (unless there are singularities 
on the curves along which we perform analytic continuation) step by step define 
holomorphic functions Xi(t), X 2 {t), a;„(t). For example, we obtain a;i(t) 

from pi{xo{t),xi{t)) = 0, then X 2 (t) from P 2 {xQ{t), xiit), X 2 (t)) = 0, etc. 

Compute the target values Xi(l), a; 2 (l), ..., Xn{l) from the given polynomial 
system, all initial values and the first target value. 

Remark 5.2. Any algorithm for this problem has to face the following funda¬ 
mental difficulty: Among all paths Xj (t) along which we perform analytic con¬ 
tinuation to define the next path Xk(t), generally only xo(t) is linear. The other 
paths xi{t), X 2 {t), ..., Xn{t) are almost always curvilinear—and unknown. We 
can at best evaluate Xi{t), X 2 {t), ..., Xn{t) at finitely many discrete points in 
time and interpolate between the sample points. However, we must make sure 
that the approximate paths we obtain by discretization remain close enough to 
the actual paths such that they yield the same result w.r.t. analytic continu¬ 
ation. In particular, we must make sure that in every step no singularities lie 
between approximate and actual path. To make things worse, this includes sin¬ 
gularities of variables that occur only in later equations, whose position in time 
may change depending on how we approximate the current step. 

One way to attack this difficulty is to eliminate cci, X 2 ^ ■ ■ ■, Xn-i from the 
polynomial system 0. e.g. using resultants. However, this approach is ex¬ 
pensive and suffers from exponential expression swell. The resulting polynomial 
equation in xq, Xn very likely has a high total degree, huge coefficients, and 
many (artificial) critical points. This means that we can in principle apply the 
method for analytic continuation of plane algebraic curves of Algorithm 3.2| but 
that in practice it will often be too expensive (see Example 7.4). If elimination 
introduces artificial critical points on the path of xp, [Algorithm 3.2| does not 
even terminate. 

Instead we pursue the following idea: 

Remark 5.3 (General scheme for homotopy continuation of triangular systems). 
We perform homotopy continuation of one equation after another, interpolating 
linearly between sample points (using a time parameter in the unit interval). 
In order to obtain sample points on the actual paths of the variables, we syn¬ 
chronize the time step. This means that we let all variables make time steps of 
the same size. We determine a step width such that analytic continuation by 


proximity is possible (as in Remark 2.2), and such that the linearly interpolated 
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paths between consecutive sample points are equivalent to the actual paths of 
the variables w.r.t. analytic continuation. To fulfil the latter requirement, the 
step width must be small enough such that there are no singularities between 
linearly interpolated paths and actual paths. We cannot foresee whether linear 
paths and actual paths enclose singularities of variables that occur only in later 
equations. We must determine whether this is the case when we later analyt¬ 
ically continue the respective variable. Should we find that we have ‘caught’ a 
singularity, we start over with a smaller step width. Unless there are singularit¬ 
ies on the actual paths of variables, there is a small neighbourhood around the 
actual paths that is free of singularities. Eventually, after finitely many reduc¬ 
tions of step size, the linear paths approximate the actual paths of the variables 
well enough such that we do not encounter singularities anymore. Then we can 
make one synchronized time step with all variables. We proceed until we reach 
time t = 1. 


6 Certified homotopy continuation of systems of 
plane algebraic curves 


In full generality, it may be very difficult to decide whether or not there are 
singularities between linearly interpolated paths and actual paths. (Among 
other things, we may want to ensure that the {k — l)-dimensional discriminant 
locus of variable Xk w.r.t. equation pk{xQ, a;i,..., Xk) = 0 does not intersect the 
polydisc around the last sample point with radii lengths of the linear paths.) 
Therefore, we restrict ourselves to systems of plane algebraic curves, a special 
case of IProblem 5.11 


Problem 6.1. Consider a system of bivariate polynomials, without loss of gen¬ 
erality 

Pi{xf),xi) = 0 , 

P2{xi,X2) = 0 , 

( 11 ) 


PniXn—l^Xji) — 0 . 

Let xo(0), a;i(0), ..., a;„(0) be initial values that satisfy the system of equations 
and let a;o(l) be a target value. We define function a;o(t) as a parameterization 
of the segment between a;o(0) and a;o(l). 


xo{t) = (1 - t)xo{0) + teo(l). 


By analytic continuation w.r.t. t € [0,1] we can (unless there are singularities 
on the curves along which we perform analytic continuation) step by step define 
holomorphic functions Xi(t), X 2 (t), ■■■, Xn(t)- For example, we obtain Xi(t) 
from pi(a;o(t), xi(t)) =0, then X 2 (t) from p 2 (xi(t), X 2 (t)) = 0, etc. 

Compute the target values a;i(l), ^ 2 ( 1 ), • ■ ■, x„(l) from the given polynomial 
system, all initial values and the first target value. 

Before we describe an algorithm for [Problem 6.11 we need the following lemma. 
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Lemma 6.2. Let C: f{x,y) = 0, xi G C be defined as in Theorem 2.1 Let 
e > 0. Suppose that we have determined S > 0 by the epsilon-delta bound of 
\ Theorem 2.1\ such that 


lyfixi) -yj{x 2 )\ < e 


for all holomorphic functions yj{x) that satisfy f{x,yj{x)) = 0 in a neighbour¬ 
hood of xi and for all X 2 with \xi — X 2 I < S. 

Then for all X 2 with S' = |xi — a; 2 | < S, 


s' = — ■ s < e 
d 


satisfies 


IVjixi) - yj{x2)\ < s'. 

This means that we can find a better estimate for the range of yj{x) w.r.t. an 
actual feasible movement of x from Xi to X 2 . 

Proof. Under the assumptions of|Lemma 6.2[ 


^ yjiSx + Xi)-yj{xi) 


is a holomorphic function from the open unit disk to the open unit disk. By the 
maximum modulus principle, we know that there exists a point on the boundary 
of the disk of radius 6' around xi where \yj(x) — yj{xi) \ is greater or equal than 
at any point x with \x — Xi\ < S'. Hence, there exists a point on the boundary 
of the disk of radius y around the origin where |/{a;)| is greater or equal than 
at any point x with |a;| < y. Schwarz lemma states that 


\f{x)\ < kl 


for all X in the open unit disk. Therefore 


e 


max 

tG[0.1] 



< max 

S' 

J 

tG[0.1] 

S 


T’ 


and thus 



for all X 2 with \xi — X 2 \ < S' < S. □ 


Remark 6.3. Alternatively, if we plug in 5 = 5' and e = e' into Q and 
s', we obtain 



MS' \ 
p{p-S')) 


< £, 


solve for 


with M, p, y as in the proof of Theorem 2.1[ This yields another better 
for the range of yj{x) w.r.t. an actual feasible movement of x from xi 


estimate 
to X 2 . 


Algorithm 6.4. Consider the system of bivariate polynomials of |Problem 6.1 
with initial values a:o(0), a;i(0), ..., a::„(0) and a target value a:o(l). 


1. Define xo{t) = (1 — t)a;o(0) + teo(l)- 
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2. Let T = 1. 

3. Let e'o = |xo(0) - a;o(r)|. 

4. For all fc = 1, 2,..., n: 

(a) Let Sk be half the minimum distance between the Xk with 

Pkixk-iiO),Xk) = 0 . 


(b) 

(c) 

(d) 

(e) 

(f) 


Compute Sk according to the epsilon-delta bound of Theorem 2.1 for 
f{x,y) = pkixk-i,Xk), xi = Xk-i{0) and £ = Efc. 

If Sk < £'k-i then let T = T/2 and go to 3. 

Let Xk(T) be the Xk with pk{xk-i(T),Xk) = 0 closest to a;fc(0). 


Let = |xfc_i(0) - a;fc_i(T)|. 

Let e'j, = ((5^ + €)/Sk • £/c with e > 0. 


5. If T = 1 then output x\{T),X 2 {T),... ,Xn{T) and stop. 

6 . Let a;o(0) = xo{T), a;i(0) = xi{T), ..., a:„(0) = Xn{T) and go to 1. 


Theorem 6.5. If the target values a;i(l), X 2 {I), ■■■, a;„(l) of Problem 6.1 are 
well-defined, \Algorithm 6.4\ computes them in finitely many steps. 


Proof. The first two steps of [Algorithm 6.4| are initialization steps. In step 1, 
we define a linear homotopy between initial value a;o(0) and final value xo(l) of 
xq. We first want to test whether we can perform analytic continuation of the 
system in a single time step. Therefore, in step 2, we set target time T = 1. 

Steps 3-6 form the main loop of our algorithm. They are repeated until we 
reach time T = 1, in which case step 5 terminates the algorithm. 

In step 3, we estimate the range of Xq as it runs from its initial position 
a;o(0) to its target position xo{T). Since Xo(t) is linear by definition (step 1), 
our estimate £q = |a;o(0) — xo{T)\ is exact. 

Step 4 is the inner loop of our algorithm, in which we try to perform analytic 
continuation equation by equation of our system. Run variable k denotes the 
index of the equation pk{xk-i,Xk) under consideration. 

In steps 4a-4b, we use the epsilon-delta bound of [Theorem 2.1| and [Re-[ 
mark 2.2 to compute a feasible step width Sk for variable Xk-i. If Xk-i moves at 


most Sk then we can perform analytic continuation of Xk w.r.t. pk{xk-i, Xk) = 0 
by selecting as Xk{T) the value of Xk with pk{xk-i{T),Xk) = 0 closest to 3:^(0). 

Hence, in step 4c, we test whether feasible step Sk is smaller than an upper 
bound of the range of Xk-i as it runs from a:fc_i(0) to a;fc_i(T). 

If Sk < we cannot be sure that there are no singularities between the 

actual path of Xk-i and the interpolated path, i.e. the segment from Xfc_i(0) to 
Xk-i{T). Our attempt to reach target time T in one step has failed. Therefore, 
we halve target time T and go back to step 3. 

Otherwise, if Sk > e'k-n tbe epsilon-delta bound of ' 


Theorem 2.1 


guarantees 

that actual path and interpolated path of Xk-i are equivalent w.r.t. analytic 
continuation of Xk. Then in step 4d, we determine target value Xk(T). By 
construction, Xk{T) is a point on the actual path of Xfc. 
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In steps 4e-4f, we use |Lemma 6.2| to compute an upper bound for the range 
of Xk as it runs from Xfc(O) to Xk{T). The computation is independent of whether 
Xk-i runs along actual or interpolated path. The bound ej, holds for both paths, 
particularly also for analytic continuation of Xk along the actual path of Xk-i- 
We then proceed with analytic continuation of the next variable, if any. 
When we leave the inner loop (step 4), we obtain valid positions for xi, a;2 ,..., 
at target time T. liT = 1, we output the solution and stop (step 5). Otherwise, 
we use xo(T),xi(r),... ,x„(T) as a valid initial position from which we again 
try to reach target time T = 1 (step 6). 

By the assumption that Xi(l), X2 (l),..., x„(l) are well-dehned, there are 
only finitely many singularities in a neighbourhood of the actual paths of Xq, 
xi, ..., x„. The algorithm terminates after finitely many steps as eventually the 
interpolated paths of xi, X2, ..., x„ approximate the actual paths well enough 
such that we do not encounter singularities anymore. □ 


7 Comparison with other approaches 


Let us discuss more examples, which allow us to compare the performance of 
our algorithm with that of other approaches. (The number of steps needed 
by [Algorithm 3.2| and [Algorithm 6.4| stated below relate to an experimental 
implementation in Haskell that is available in the supplementary material for 
this article.) 

Example 7.1 ( Hauenstein et al. (2014[ Section 7.1)). Consider the Newton 
homotopy 

H (x, t) = /(x) + vt 

where /(x) = x^ — 1 — to and v = m for various values of to > —1. The goal is 
to analytically continue x as t moves from 1 to 0. 

In [Table 7.^ and [Table 7.3 we compare the performance of [Algorithm 3.2 
with that of the algorithms of Beltran and Leykin (2013) and Hauenstein et al. 
[(2014[), for various values of to. The data for the latter algorithms is quoted 


from ([Hauenstein et ah, 2014 Table 1 and Table 2). 

Both [Beltran and Leykin (2013| ) and Hauenstein et al. (2014| | present algorithms 
designed for certified homotopy continuation of arbitrary polynomial systems 
whereas [Algorithm 3.2[ can only deal with plane algebraic curves. However, the 
example indicates that in the univariate case [Algorithm 3.2| may perform much 
better than those more general algorithms, which do not exploit the special 
structnre of the univariate case. 


Furthermore, let us elaborate on Remark 5.2[ The following example shows that 
it may be better to apply [Algorithm 6.4 to a system of plane algebraic curves 
than to eliminate variables and apply Algorithm 3.2[ to the resultant. 

Example 7.4. Consider the system of bivariate polynomials 
Pi(xo, xi) = —4 + 2xo + xi + 2xoXi + x^ = 0, 

P 2 iXi,X 2 ) = xf + X 2 = 0, 
with initial values 


( 12 ) 


xo(0)=0, xi(0) = 


-1 - Vi7 


2:2(0) = 


'-9- 
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Number of steps 
of |Algorithm 3.2] 


Number of a priori 
steps of|Beltran 
Leykin (2013) 


Number of a pos¬ 
teriori certified in¬ 
tervals of IHauen-l 
stein et al. (20141 


10 

9 

184 

51 

20 

12 

217 

67 

30 

14 

237 

78 

40 

16 

250 

82 

50 

17 

260 

88 

60 

18 

269 

92 

70 

19 

276 

96 

80 

20 

282 

99 

90 

21 

288 

103 

100 

21 

292 

105 

1000 

41 

395 

162 

2000 

49 

426 

180 

3000 

54 

446 

191 

4000 

58 

457 

197 

5000 

62 

468 

204 

10000 

73 

499 

220 

20000 

87 

530 

238 

30000 

96 

547 

250 


Table 7.2: Performance of Algorithm 3.2 in comparison with the algorithms 
of Beltran and Leykin (2013) and Hauenstein et al. (2014), for various values 
of m. The data in the last two columns is quoted from |Hauenstein et al., 


2014 


Table 1). 


Number of steps 
of [Algorithm 3.2 


Number of a priori 
steps ofjBeltran and 
Leykin (2013) 


Number of a pos¬ 
teriori certified in¬ 
tervals of IHauen-l 
stein et al. (2014) 


1 

5 

176 

64 

2 

9 

287 

68 

3 

14 

390 

70 

4 

18 

492 

71 

5 

22 

593 

71 

6 

27 

695 

71 

7 

31 

798 

71 

8 

36 

901 

71 

9 

40 

1003 

71 

10 

44 

1108 

71 


Table 7.3: Performance of 

Algorithm 3.2 

in comparison with the algorithms 

of 

Beltran and Leykin (2013 

) and 

Hauenstein et al. (2014), for various values of 


et ah, 2014 Table 2). 
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and target value a:o(l) = 1. The xi-resultant of pi(a;o,xi) and ^ 2 ( 2 ^ 1 , 2 : 2 ) is 

q{xo, X 2 ) = 16 — 16x0 + 4xo + ^^2 + 4x0^2 + xf = 0. (13) 


Let us compare the performance of Algorithm 6.4 for (12) with the performance 


of Algorithm 3.2 for (131 as xq moves linearly (in the unit interval) from 0 to 1. 
We find that [Algorithm 6 "^ subdivides once, i.e. it needs two steps. In contrast, 
Algorithm 3.2 needs six steps. One possible explanation is that Xq = — 5 is 


a singularity of (13) but not of (12). For xq = — i, (13) has three zeros of 


multiplicity two, whereas (12) has six simple roots. Each zero of multiplicity 


two of (13) corresponds to two simple zeros of (12) with differing signs of Xi. 

Generally, elimination introduces artificial singularities. Due to an artifi¬ 
cial singularity it can even happen that we cannot analytically continue the 
resultant: For example, the xi-resultant of 

Pi(xo, xi) = — 4 -I- 2xo -I- X 2 — 2xoXi -I- x^ = 0, 

P 2 (a:i,X 2 ) = Xi -I- X 2 = 0, 


has an artificial singularity at Xq = |. In this case. Algorithm 3.2 
terminate whereas [Algorithm 6 . 4| produces the desired result. 


does not 


8 Conclusion 


From an epsilon-delta bound for plane algebraic curves (Theorem 2.1), we have 


derived algorithms for certified homotopy continuation of plane algebraic curves 
(Algorithm 3.2) and systems of plane algebraic curves (Algorithm 6.4). Our 


certificate is rigorous for exact real arithmetic. For floating point arithmetic, 
[Theorem 2.1| can be considered a soft certificate. Several examples demonstrate 
the practicability of our approach. 

A generalization of [Algorithm 6.4' to arbitrary systems of polynomials might 
be an interesting challenge for further research. 


Acknowledgements 

I thank Ulrich Bauer, Tim Hoffmann, Jurgen Richter-Gebert, Katharina Schaar, 
and Martin von Gagern for their support in preparing this article. 


Funding 

This research was supported by DFG Collaborative Research Center TRR 109, 
“Discretization in Geometry and Dynamics”. 


References 

Ahlfors, Lars Valerian. 1979. Complex Analysis, 3rd ed., McGraw-Hill, Singapore. 

Allgower, Eugene L. and Kurt Georg. 1990. Numerical Continuation Methods: An Introduc¬ 
tion, Springer Series in Computational Mathematics, vol. 13, Springer, Berlin. 

Beltran, Carlos and Anton Leykin. 2012. Certified Numerical Homotopy Tracking, Experi¬ 
mental Mathematics 21, no. 1, 69-83, DOT 10.1080/10586458.2011.606184. 


19 
































_2013. Robust certified numerical homotopy tracking, Foundations of Computational 

Mathematics 13, no. 2, 253-295, DOT 10.1007/sl0208-013-9143-2. 

Fujiwara, Matsusaburo. 1916. Uber die obere Schranke des absoluten Betrages der Wurzeln 
einer algebraischen Gleichung, Tohoku Mathematical Journal, First Series 10, 167-171. 

Hauenstein, Jonathan D., Ian Haywood, and Alan C. Liddell Jr. 2014. An a posteriori cer¬ 
tification algorithm for Newton homotopies, ISSAC ’14 (Kobe, Japan, 2014), Proceedings 
of the 39th International Symposium on Symbolic and Algebraic Computation, ACM, New 
York, NY, USA, pp. 248-255, DOT 10.1145/2608628.2608651, (to appear in print). 

Hauenstein, Jonathan D. and Frank Sottile. 2012. Algorithm 921: alphaCertified: Certifying 
Solutions to Polynomial Systems, ACM Transactions on Mathematical Software 38, no. 4, 
28:1-28:20, DOT 10.1145/2331130.2331136. 

Hoffmann, Tim. 2009. Discrete Differential Geometry of Curves and Surfaces, MI Lecture 
Notes, vol. 18, Faculty of Mathematics, Kyushu University, Japan. 

Kortenkamp, Ulrich. 1999. Foundations of Dynamic Geometry, dissertation, ETH Zurich, 
Zurich. 

Kortenkamp, Ulrich and Jurgen Richter-Gebert. 2001. Grundlagen dynamischer Geomet- 
rie, Zeichnung - Figur - Zugfigur: Mathematische und didaktische Aspekte dynamischer 
Geometrie-Software (H.-J. Elschenbroich, Th. Gawlick, and H.-W. Henn, eds.), Franzbecker, 
Hildesheim, pp. 123-144. 

_ 2002. Complexity issues in dynamic geometry, Festschrift in the honor of Stephen 

Smale’s 70th birthday (M. Rojas and Felipe Gucker, eds.), World Scientific, pp. 355-404. 

_ 2006. Cinderella: The interactive geometry software, http://www.cinderella.de 

Morgan, Alexander and Andrew Sommese. 1987. A Homotopy for Solving General Polyno¬ 
mial Systems That Respects m-Homogeneous Structures, Applied Mathematics and Compu¬ 
tation 24, 101-113. 

Smale, Steve. 1986. Newton’s Method Estimates from Data at One Point, The Merging of 
Disciplines: New Directions in Pure, Applied, and Computational Mathematics (Richard E. 
Ewing, Kenneth 1. Gross, and Clyde F. Martin, eds.), Springer, New York, pp. 185-196, DOI 
10.1007/978-l-4612-4984-9_13, (to appear in print). 

Sommese, Andrew J. and Charles W. Wampler H. 2005. The Numerical Solution of Systems 
of Polynomials Arising in Engineering and Science, World Scientific, Singapore. 


20 





